COMMENT ⊗   VALID 00007 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	TITLE FITS    -  A LEAST SQUARES CUBIC FIT  -  NOVEMBER 1972.
C00004 00003	ACCUMULATE PRODUCTS OF THE COORDINATES OF THE DATA POINTS.
C00006 00004	Solve four equations in four unknowns by Gaussian Elimination.
C00008 00005	Find column-2 pivot for 3 by 3 matriX
C00009 00006	Find column-3 pivot for 2 by 2 matriX
C00010 00007	Back Substitution.
C00011 ENDMK
C⊗;
TITLE FITS    -  A LEAST SQUARES CUBIC FIT  -  NOVEMBER 1972.

COMMENT/

   Y = (((B0*X + B1)*X + B2)*X + B3).

	  SY = B3*N   +  B2*SX +  B1*X2  +  B0*X3.
	  XY = B3*SX  +  B2*X2 +  B1*X3  +  B0*X4.
	 X2Y = B3*X2  +  B2*X3 +  B1*X4  +  B0*X5.
	 X3Y = B3*X3  +  B2*X4 +  B1*X5  +  B0*X6.
/

;COEFFICIENTS OF THE FOUR EQUATIONS.

	ROW1:	BLOCK 6
	ROW2:	BLOCK 6
	ROW3:	BLOCK 6
	ROW4:	BLOCK 6

;Use Accumulators for accumulating.
	ACCUMULATORS{X,Y,SX,SY,X2,XY,X3,X2Y,X4,X3Y,X5,X6}
		I←1

;CUBFIT(A,B,N)
; - halfword array of N points (SX,y) given in A.
; - returns four coefficients in real array B.

SUBR(CUBFIT)

	dac 12,AC12
	lac I,ARG3
	lacn ARG1↔SOS
	hrlm 0,I		;loop count and pointer.
	lac[xwd 4,5]↔setz 4,↔blt 15	;clear accumulators.
;ACCUMULATE PRODUCTS OF THE COORDINATES OF THE DATA POINTS.

L1:	NIP X,(I)	↔	FSC X,233-12
	NAP Y,(I)	↔	FSC Y,233-12

	FADR SX,X
	FADR SY,Y
	FMPR Y,X	↔	FADR XY,Y
	FMPR Y,X	↔	FADR X2Y,Y
	FMPR Y,X	↔	FADR X3Y,Y

	LAC X
	FMPR X↔FADR X2,
	FMPR X↔FADR X3,
	FMPR X↔FADR X4,
	FMPR X↔FADR X5,
	FMPR X↔FADR X6,

L2:	AOBJN I,L1

;SETUP THE COEFFICIENTS.

	LAC 0,ARG1↔AOS↔FSC 0,233
	MOVEI 1,ROW1
	MOVEI 2,ROW2
	DAC  0,(1)1↔DAC SX,(1)2↔DAC X2,(1)3↔DAC X3,(1)4↔DAC SY,(1)5
	DAC SX,(2)1↔DAC X2,(2)2↔DAC X3,(2)3↔DAC X4,(2)4↔DAC XY,(2)5
	MOVEI 3,ROW3
	MOVEI 4,ROW4
	DAC X2,(3)1↔DAC X3,(3)2↔DAC X4,(3)3↔DAC X5,(3)4↔DAC X2Y,(3)5
	DAC X3,(4)1↔DAC X4,(4)2↔DAC X5,(4)3↔DAC X6,(4)4↔DAC X3Y,(4)5
;Solve four equations in four unknowns by Gaussian Elimination.

;TRIANGULARIZATION.
	;Accumulators 1,2,3,4 are row pointers.
	R←5	;Reciprocal of the Pivot.
	M←6	;Multiplier of the row.
	SETZM PARITY#

;Find column-1 pivot for 4 by 4 matrix.

T1:	lacm (1)1↔lacm M,(2)1↔caml M↔jrst .+3↔exch 1,2↔setcmm PARITY
	lacm (1)1↔lacm M,(3)1↔caml M↔jrst .+3↔exch 1,3↔setcmm PARITY
	lacm (1)1↔lacm M,(4)1↔caml M↔jrst .+3↔exch 1,4↔setcmm PARITY

;Reciprocal of the Pivot.

	movsi R,(1.0)↔fdvr R,(1)1

;Zero column-1 elements of rows 2, 3, & 4.

	lac M,(2)1↔fmpr M,R↔	 SETZM (2)1
		lacn M↔fmpr (1)2↔fadrm (2)2
		lacn M↔fmpr (1)3↔fadrm (2)3
		lacn M↔fmpr (1)4↔fadrm (2)4
		lacn M↔fmpr (1)5↔fadrm (2)5

	lac M,(3)1↔fmpr M,R↔	 SETZM (3)1
		lacn M↔fmpr (1)2↔fadrm (3)2
		lacn M↔fmpr (1)3↔fadrm (3)3
		lacn M↔fmpr (1)4↔fadrm (3)4
		lacn M↔fmpr (1)5↔fadrm (3)5

	lac M,(4)1↔fmpr M,R↔	 SETZM (4)1
		lacn M↔fmpr (1)2↔fadrm (4)2
		lacn M↔fmpr (1)3↔fadrm (4)3
		lacn M↔fmpr (1)4↔fadrm (4)4
		lacn M↔fmpr (1)5↔fadrm (4)5

;Normalize First Row.

	FMPRM R,(1)1
	fmprm R,(1)2
	fmprm R,(1)3
	fmprm R,(1)4
	fmprm R,(1)5
;Find column-2 pivot for 3 by 3 matriX
T2:	lacm(2)2↔lacm M,(3)2↔caml M↔jrst .+3↔exch 2,3↔setcmm PARITY
	lacm(2)2↔lacm M,(4)2↔caml M↔jrst .+3↔exch 2,4↔setcmm PARITY

;Reciprocal of the pivot.

	movsi R,(1.0)↔fdvr R,(2)2

;Zero column-2 elements of rows 3 & 4.

	lac M,(3)2↔fmpr M,R	↔SETZM (3)2
		lacn M↔fmpr (2)3↔fadrm (3)3
		lacn M↔fmpr (2)4↔fadrm (3)4
		lacn M↔fmpr (2)5↔fadrm (3)5

	lac M,(4)2↔fmpr M,R	↔SETZM (4)2
		lacn M↔fmpr (2)3↔fadrm (4)3
		lacn M↔fmpr (2)4↔fadrm (4)4
		lacn M↔fmpr (2)5↔fadrm (4)5

;Normalize Second Row.

	FMPRM R,(2)2
	fmprm R,(2)3
	fmprm R,(2)4
	fmprm R,(2)5
;Find column-3 pivot for 2 by 2 matriX

T3:	lacm(3)3↔lacm M,(4)3↔caml M↔jrst .+3↔exch 3,4↔setcmm PARITY

;Reciprocal of the pivot.

	movsi R,(1.0)↔fdvr R,(3)3

;Zero column-3 element of row 4.

	lac M,(4)3↔fmpr M,R	↔SETZM (4)3
		lacn M↔fmpr (3)4↔fadrm (4)4
		lacn M↔fmpr (3)5↔fadrm (4)5

;Normalize Row-3.

	FMPRM R,(3)3
	fmprm R,(3)4
	fmprm R,(3)5

;Find column-4 pivot for 1 by 1 matriX  (trivail).
;Reciprocal of the pivot.

T4:	movsi R,(1.0)↔fdvr R,(4)4

;Normalize Row-4.

	FMPRM R,(4)4
	fmprm R,(4)5
;Back Substitution.
	
T5:	lacn (4)5↔fmpr (3)4↔fadrm (3)5
	lacn (4)5↔fmpr (2)4↔fadrm (2)5
	lacn (4)5↔fmpr (1)4↔fadrm (1)5
	lacn (3)5↔fmpr (2)3↔fadrm (2)5
	lacn (3)5↔fmpr (1)3↔fadrm (1)5
	lacn (2)5↔fmpr (1)2↔fadrm (1)5

;Move results to Vector B.
	lac  5,ARG2	;B vector pointer.
	skipe PARITY↔jrst T7

T6:	lac(1)5↔dac(5)3
	lac(2)5↔dac(5)2
	lac(3)5↔dac(5)1
	lac(4)5↔dac(5)0
	lac 12,AC12↔pop3j
	
T7:	lacn(1)5↔dac(5)3
	lacn(2)5↔dac(5)2
	lacn(3)5↔dac(5)1
	lacn(4)5↔dac(5)0
	lac 12,AC12↔pop3j
END